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A good representation of mesoscopic fluids is required to combine with molecular simulations at 
larger length and time scales (De Fabritiis et. al, Phys. Rev. Lett. 97, 134501 (2006)). However, 
accurate computational models of the hydrodynamics of nanoscale molecular assemblies are lacking, 
at least in part because of the stochastic character of the underlying fluctuating hydrodynamic 
equations. Here we derive a flnite volume discretization of the compressible isothermal fluctuating 
hydrodynamic equations over a regular grid in the Eulerian reference system. We apply it to fluids 
such as argon at arbitrary densities and water under ambient conditions. To that end, molecular 
dynamics simulations are used to derive the required fluid properties. The equilibrium state of the 
model is shown to be thermodynamically consistent and correctly reproduces linear hydrodynamics 
including relaxation of sound and shear modes. We also consider non-equilibrium states involving 
diffusion and convection in cavities with no-slip boundary conditions. 

PACS numbers: 47.61.-k,47.11.Mn 



I. INTRODUCTION 

^From a continuum perspective the fundamental equa- 
tions underlying hydrodynamics at the mesoscale are the 
well known fluctuating hydrodynamics (FH) equations 
The FH equations are stochastic partial differential 
equations that reduce to the Navier-Stokes equations in 
the limit of large volumes. In fact, at scales of nano to mi- 
cro meters, thermal fluctuations cannot be neglected, but 
must be incorporated as random terms in the momentum 
and energy equations of hydrodynamics because they are 
responsible for the mechanical and thermal energy pro- 
cesses underlying Brownian motion. The hydrodynamics 
of nanoscopic quantities of liquids, in particular for water 
models used in molecular dynamics simulations such as 
TIP3P 0], is relevant in many biological and technolog- 
ical applications. 

Novel multiscale modelling techniques via domain de- 
composition (particle-continuum hybrid approach) re- 
quire a very accurate description of thermodynamics and 
hydrodynamics at the mesoscale level [1, In these 
methods a large part of the system is resolved with a con- 
tinuum model (CFD) and a smaller part using full-atom 
molecular dynamics (MD). An exact match between the 
local thermodynamic and hydrodynamic properties of 
the continuum and the molecular system is required for 
such a scheme to work properly. This matching enables 
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a seamless coupling such that the behaviour of molecu- 
lar and hybrid simulations are indistinguishable 4]. De- 
pending on the process or regime being considered fluctu- 
ations may play an important role For these reasons, 
the fluctuating hydrodynamics model described here has 
been used in the first hybrid MD-FH model for water 
which includes mass and momentum fluctuations as well 
as for propagation of sound waves across an hybrid in- 
terface [1,11] because of the need to have a fine control 
over the characteristics of the fluid to much the molec- 
ular description (shear and bulk viscosities, equation of 
state and thermodynamic fluctuations). 

An accurate code for fluctuating hydrodynamics, pos- 
sibly interfaced with molecular dynamics, would be a 
useful tool for nanoscale computational fluid dynamics 
(CFD) simulation, including inter alia microfluidic de- 
vices These devices are essentially hydraulic mi- 
cro machines which are able to process nano liters of 
reagents. These volumes are too large to be simulated 
by molecular dynamics, while, on the other hand, a stan- 
dard CFD code cannot handle fluctuations at all. 

A general purpose FH solver could also be used to 
provide an implicit hydrodynamic solvent for solute par- 
ticles (polymers, colloids, etc.). Solvent molecules of- 
ten comprise the computationally most expensive part of 
any molecular or coarse-grained simulation, but in some 
instances the solvent could be approximated by an "im- 
plicit description" , retaining only the hydrodynamic con- 
tribution to the solute dynamics. A possible approach 
was flrst illustrated in [sj . A solver of FH for the hydro- 
dynamic description of these hybrid models could be em- 
ployed to study the effect on polymer collapse of changing 
solvent characteristics @ by tuning the characteristics of 
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the fluids such as their viscosities. 

In recent years, many computational models have 
been devised which provide a discrete representation of 
FH [lol [Til e. g. dissipative particle dynamics 
(DPD) Imllland the lattice Boltzmann method (LBM) 
[Til . [Ta . TTgI. IITL [T^ later extended to include thermody- 
namic fluctuations [T^, [Ml . Following a continuum 
approach, the fundamental equations of fluctuating hy- 
drodynamics can also be resolved directly via flnite differ- 
ences or finite volume schemes [22!, . Most of the pre- 
vious schemes proposed for solving fluctuating hydrody- 
namics considered only the gas phase [H, [l^] and have fo- 
cused on the linear regime where a closed set of equations 
for the fluctuating quantities or for their mutual spatial 
correlations can be derived The main difficulty in 

devising a discretization of the FH equations is that the 
precise form of the required fluctuation-dissipation rela- 
tions depends on the discretization scheme and, in gen- 
eral, it does not coincide with the fluctuation-dissipation 
relations of the continuum description. This statement 
applies for Lagrangian FH models f23', '25', '26*1 , for fluctu- 
ating lattice Boltzmann models JJ,, 2ih 21\ and for the 
Eulerian FH description. Moreover the resulting equa- 
tions are stochastic in nature which adds extra compli- 
cations to the integration methods i28|, i2^] . 

For these reasons, a simple and easy-to-implement gen- 
eral purpose solver of the compressible FH equations al- 
lowing fluid specificity and non-linear hydrodynamic cou- 
pling is not readily available. It is the purpose of this 
paper to address this lacuna. We present a finite volume 
discretization of the compressible isothermal fluctuating 
hydrodynamic equations, based on an Eulerian descrip- 
tion on a regular grid. The model provides a thermo- 
dynamically consistent coarse-grained representation of 
nano liter portions of real fluids ranging from gases (ar- 
gon) to liquids (water). 

The outline of the paper is the following: The fluctuat- 
ing hydrodynamic Eulerian solver is described in section 
im In section Hill we study the numerical accuracy of the 
scheme by comparing the input values of the viscosities 
with the effective ones measured from the hydrodynamic 
solver based on the relaxation of sound and shear waves. 
In section IIVI we assess the validity of the description 
of the equilibrium state, showing that fluctuations are 
correctly generated, propagated and dissipated. To that 
end we calculate the time correlation functions of the 
different fluctuating variables of one fluid cell (density 
and velocity) and compare them with the corresponding 
grand canonical result. In section |V] we consider non- 
equilibrium states in closed systems with rigid walls, us- 
ing the no-slip boundary condition: we test them for Cou- 
ette, Poiseuille and cavity flows. Finally, we summarise 
our findings in section IVll 



II. A FINITE VOLUME DISCRETIZATION OF 
FLUCTUATING HYDRODYNAMICS 

Our proposed mesoscopic model is a finite volume dis- 
cretization of FH [2I, [201 over a regular lattice in the 
Eulerian frame of reference. In this case, we concentrate 
on the description of an isothermal compressible fiuid. 
This sort of description can be generalized straightfor- 
wardly to non-isothermal states for fluids with vanish- 
ingly small thermal expansion, for which the energy equa- 
tion decouples from the mass and momentum equations 
[23 |. Extensions to include energy flows will be consid- 
ered elsewhere. We thus require the equations of fluctu- 
ating hydrodynamics describing the conservation of mass 
and momentum, 

dtp = -dpgp, 

dtQa = -dp (^gfjVa + Ha/3 + Ha/j) , (1) 

where p{r,t) is the density field of the fluid, Va{r,t) 
is the continuous velocity fleld in the component a, 
gp{r,t) = p{r,t)vp{r,t) is the momentum field and we 
have used the repeated suffix convention for summation 
over repeated indices. Uap{r, t) and Iiap{'^^ t) are respec- 
tively the average (Navier-Stokes) and fiuctuating stress 
tensor fields. The average stress tensor is defined as 
n = (j5-|-7r)H-n, where p is the thermodynamic pressure 
given by the equation of state for the fluid, tt = —C,d~fV~f 

and liap — —11 {daVp -t- dpVa — 2D~^d^v.ySap) where 77 
and C are the shear and bulk viscosity respectively and 
D is the spatial dimensionality. 

The equation of state p — p{p, T) for Lennard- Jones 
(LJ) fluids (like argon) has been studied by several au- 
thors (see e.g., Ref. I30i1 ) , as well as the transport coeffi- 
cients of the L J fluid |3ll. Is^ . By contrast, the equation 
of state of the TIP3P model for water [2] (chosen due 
to its importance in biological applications) has received 
less attention [33| . In Appendix [B] we provide a para- 
metric study of the equations of state of the fluids consid- 
ered here (argon and TIP3P water model), performed via 
molecular dynamics (MD) simulations. From this study 
we obtain a second order polynomial fit for p = p(p, T) 
which provides the equation of state of our FH model. 
This procedure is required, for instance, to provide the 
level of accuracy for the thermodynamic pressure within 
hybrid MD simulations Q. We have also calculated the 
transport coefficients of water via non-equilibrium molec- 
ular dynamics, which are in agreement with those re- 
ported in the literature [3^ . 

The fiuctuating stress tensor Ilap (see Ref. [J ) is a ran- 
dom Gaussian matrix with zero mean and correlations 
given by 

(na/3(ri,ti)n57(r2,i2)) ^2kBTCai3js5{ti -t2)S{ri -r2), 

(2) 

where CaPjS = [viSasSfBj + Sa^SiSg + (C - J)'n)^al3Ss'i] , 

ks is the Boltzmann constant and T is the temperature. 
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Note that this spatial delta-correlated quantity, in the 
discrete limit of a small volume and small time interval, 
can be rewritten as 



long algebra is avoided; we find discrete versions for 



2kBT_ 
'AtAV 
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where AT^ is the small volume element of fluid and At is 
the time step. 

The fluctuating hydrodynamic equations ([T]) are bal- 
ance equations of the form 9t(/)(r, t) — —V • J"^ for mass 
and momentum which can be integrated by considering 
a finite volume discretization. In what follows we will 
derive a finite volume discretization of the equations of 
fluctuating hydrodynamics in the Eulerian system of ref- 
erence. We first partition the space into N space filling 
volumes Vfe (in our case a regular Cartesian lattice is 
used) with k = 1,...,N to integrate Eqs. ^ over the 
volume Vk and apply Gauss's theorem 



di 
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where e^,; is the unit vector perpendicular to the con- 
tact surface of area Aki from volume I to volume k. The 
summation is over all the I volumes that are in contact 
with fluid volume k. By defining Mj^ = Jy^p{r,t)dr 

as the mass inside a generic volume Vk, and P|, = 
Jy^ p(r, t)v(r, i)(ir its momentum, we thus build dynam- 
ical equations corresponding to discrete extensive vari- 
ables which replace Eqs.(IT]) governing the time evolution 
of the intensive continuum fields. These new equations 
are given by 

dMl = ^gki- BkiAkidt, 
I 

dPl = ^ (vfe, gfc, +Uki)- ekiAkidt + dPl (5) 
I 



where we have approximated the mass flux J^^ with 
gki = liPk + Pi)^i'Vk + vi), the velocity on the surface 
kl as Vki = ^(vfc -I- vj), the average stress tensor on the 

surface as 11^; = ^[{pi + tt;)! + 11/] and rfP^. indicates 
the momentum change due to the fluctuating part of the 
pressure tensor, all at time t. 

An essential component for the discretization of a 
stochastic mesoscopic model is the balancing of dis- 
sipative and fluctuating components, otherwise the 
fluctuation-dissipation theorem would not be satisfied. 
There are at least two ways to satisfy this condition, ei- 
ther by using the Fokker-Planck equations mathemati- 
cally equivalent to the stochastic differential equations 
(SDE [5]) and equating dissipative and diffusive terms 
weighted over the Gibbs ensemble distribution [2^ or 
by using the GENERIC formalism [H, [H. By choos- 
ing the gradient discretization provided in a lot of 
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and the fluctuating component of the momentum equa- 
tion given by 



Ajd 
2 



(7) 



where dW; is a 13 x D matrix (13 = 3 in three di- 
mensions) of independent Wiener increments satisfying 

(dWfdWf) = SkiSa-ySpsdt and dW, is a traceless 
symmetric random matrix deflned as 



tr[dWi] ^ 
D 



(8) 



The resulting set of stochastic differential equations is 
integrated using a simple stochastic Euler scheme in the 
present work. Note, however, that other more accu- 
rate stochastic integration schemes for mesoscopic mod- 
els have recently been proposed based on the Trotter ex- 
pansion in the stochastic case [1^ . Improvements to 
the solver for the spatial regular grid and for the time in- 
tegration scheme have not been considered for the present 
scheme because at these scales the Reynolds number 
is usually low and the computational limitation comes 
rather from the molecular dynamics component. How- 
ever, generalizations to unstructured grids are straight- 
forward. In particular, in a hybrid MD-CFD model finer 
cells should be located near the MD region and coarser, 
bigger cells further away from the MD domain. 



III. ACCURACY OF THE SCHEME 

To assess the accuracy of the numerical scheme, we 
measure the effective viscosities and sound velocity com- 
puted from the hydrodynamic solver and compare them 
with the input values. The transport coefficients are mea- 
sured from the relaxation of transversal and longitudinal 
waves in the deterministic limit. We give more details 
in Appendix |^ In the following tests we consider ar- 
gon at temperature T = 300 K and mass equilibrium 
density pe = 0.6 g/mol/A^ and TIP3P water. The corre- 
sponding values of the dynamic shear and bulk viscosity 
for argon and water are shown in table [H along with the 
isothermal sound velocity, c\ = {dP/dp)T and sound ab- 
sorption coefficient Ft- We consider a periodic 3D cubic 
domain of length L; thus, the permitted wavelengths are 
kn = 27rn/L, and we excite the longest wavelength of the 
system, of wave- number 27r/L. 
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TABLE L Some properties of argon and water at T = 300 K 
at the mass densities p considered; m is the molecular mass. 
Note that length, time and mass units are A, ps and g/mol 
respectively. The properties displayed are shear viscosity (77), 
bulk viscosity [C,), isothermal sound speed (ct) and isothermal 
sound absorption (Ft) in corresponding units. 
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water 


18.015 


0.632 


53.71 
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14.75 


157.17 
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FIG. 1: Velocity field function of z (in A) for 

this decaying transversal wave. Diamonds correspond to 
simulation results and the continuum lines correspond to 
the theoretical profile at snapshots corresponding to times 
i = 0, 50, 100, • • • , 750 ps. At time t = ps the amplitude is 
maximum while at t = 750 ps the amplitud is mininum. 



A. Transversal wave 

The wave- vector of a transversal wave is perpendicular 
to its velocity, k_Lv. Consider an initial perturbative 
velocity v = (uq sin(fc2:), 0, 0), with k = (0,0, /c). For 
small perturbations around equilibrium, the linearized 
solution for the momentum density field in Fourier space 
given in Eqs. (|A9|) is 



g(fc,i) = exp{-iykh}g{k,0), 



(9) 



while in real space, the time dependent velocity field 
reads 



Vx{t) = vosin{kz) exp{—iyk^t}; 
in addition Vy{t) — Vz{t) — and p(t) ~ Pf.. 



(10) 



In Fig[T]we plot some snapshots of the velocity field for 
a deterministic simulation. This case represent a three 
dimensional simulation box of 200 x 200 x 200 (10 x 
lOx 10 cells) which corresponds to a spatial resolution 5 = 
20 A. The applied initial velocity amplitude is wo — 2.04 
A/ps, the argon input viscosity at 300 K is 77 = 5.4744 
and the mean mass density is p = 0.6 g/mol/A'^. The 
theoretical agreement with expression (jlOp is remarkable. 
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FIG. 2: Relative error for the shear viscosity, sound absorp- 
tion coefficient and sound velocity as a function of the spatial 
resolution 5. For any fluid parameter (e.g. u) the relative 
error is defined as = \vnum — v\lv, where v is the in- 
put value and Vnum that measured from the relaxation of the 
corresponding hydrodynamic mode. 



For a closer inspection of the numerical accuracy of the 
scheme we compared the effective (or numerical) shear 
viscosity Vnum with the input value v. The value of 
Vnum was measured by fitting the decay of the Fourier 
component of the transversal momentum to a simple 
exponential function. The relative error in viscosity 
Ey = \vnum — is shown in Figl2] against the spatial 
resolution (5, given by the distance between contiguous 
cells. The trend obtained is Ey oc 6^ '^'^ ^ showing that our 
spatial discretization method is of second order. As an 
example, a continuum cell size of (5 = 20 A wiU give an 
error in the viscosity around 12% while (5 = 15 A will re- 
duce the relative error to 5% for the pertrubation applied 
here. 



B. Longitudinal wave 

If the equilibrium state of the fluid is initially per- 
turbed with a momentum field go = (50 sin /ex, 0, 0), a 
longitudinal sound wave is created with a wave-vector 
parallel to the momentum perturbation fc — (27r/L, 0, 0). 
As shown in Appendix [Al two travelling sound modes 
propagate at the sound velocity ct, creating a standing 
wave in the periodic domain. According to Eqs. (|A9[) 
the Fourier components of the deviation from the equi- 
librium state for the density and velocity of the sound 
modes are given by 



6p{k,t) oc —iexp{—TTk^t}sm{cTkt) 



k ■ go 



6g{k,t) oc exp{—TTk^t} cos{cTkt)k 



CT 

k ■ go 

CT 



[11) 
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FIG. 3: The mass density p as a function of x (in A) for a de- 
caying longitudinal wave. Diamonds correspond to simulation 
results and the continuum lines correspond to the theoretical 
profile given by Eq. (|12p. The snapshots correspond to the 
times t = Q, 25, 50, ■ ■ • , 750 ps. The agreement is remarkable. 
At time t — ps the amplitude is maximum while at t = 750 
ps the amplitud is mininum. 



while in real space the density and velocity fields evolve 
like 

p{t) = Pe + — cos(fcx) sin(cTfct) exp{— r^fc^i}, 

CT 

Vx{t) = — sm{kx) cos{cTkt) expl—TTk^t} (12) 

Pe 

and the other velocity field components remain at rest. 
Note that Pe denotes the equilibrium density. 

In Figl3] we plot the density field as a function of x 
for a deterministic three dimensional simulation of size 
200 X 200 X 200 (10 x 10 x 10 cells). In order to keep 
the system in the linear regime, the applied initial veloc- 
ity amplitude is set to vq = 0.0204 A/ps and the argon 
input mean mass density is p = 0.6 g/mol/A'^. The best 
fit to the absorption coefhcient and sound velocity are 
Tt = 7.24 and ct = 5.42, providing relative errors of 
4.8% and 3.6% respectively. Figure [2] also shows the rel- 
ative error in the sound absorption coefficient and sound 
velocity versus the spatial resolution. Note that the rel- 
ative error in both quantities also decays roughly like 6^ 
which corroborates the second order spatial resolution of 
the scheme again. 



IV. THE EQUILIBRIUM STATE 

Within a fluid volume, stress fluctuations arise due to 
forces involved in the random sequence of molecular col- 
lisions. This fluctuating force generates momentum fluc- 
tuations. The amplitudes of mass and momentum fluc- 
tuations are determined by thermodynamic constraints. 
Each fluid cell is an open system with constant volume 
and temperature, in other words, belonging to the grand 
canonical ensemble. The variance of momentum and 




V (A/ps) 

FIG. 4: The equilibrium distribution of the x component of 
the velocity at one fluid cell of volume 37.5 nm^ in a simula- 
tion of argon (circles) at T = 300 K and p = 0.6 g/mol/A^ 
compared with the theoretical normal distribution (contin- 
uum line). In this simulation At = 20 fs and the best normal 
fit to the numerical distribution yields T„um = 296.28 K, that 
is a relative error of around 1.2%. 



mass of a cell "c" are [i| 

Var[P, -e] = pVcksT, 

Var[M,] = (13) 

where e is a unit surface vector, the mean mass of the cell 
is pVc and ct is the isothermal sound velocity. On the 
other hand, these spontaneous mass and momentum fluc- 
tuations are transported through the fluid and dissipated 
following the same mechanism underlying the hydrody- 
namic modes explained in Appendix [Al i.e., either via 
shear or sound modes. 



A. Amplitude of fluctuations 

In this section we consider the equilibrium state of ar- 
gon at different densities (from gas to liquid) and wa- 
ter (TIP3P model, see Appendix [B]) in order to illus- 
trate that the fluctuations are generated, transported 
and dissipated in a thermodynamically and hydrodynam- 
ically consistent way. First, we confirm that the ampli- 
tudes of mass and momentum fluctuations are consis- 
tent with thermodynamic relations plj]. In the numeri- 
cal scheme the amplitude of fluctuations is determined, 
by construction, via the fluctuation-dissipation theorem 
j3q . Figure 2] presents a typical distribution of one ve- 
locity component and compares it with the theoretical 
Maxwellian distribution. As usual, a temperature can 
be extracted from the variance of the velocity distribu- 
tion. This "numerical" temperature will be labelled as 
Tnum = X)q Var[i;Q]/(3A:B). 
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FIG. 5: Relative error in the temperature of the scheme, Et ~ 
\Tnum - T\/T, with T„um = Ec VsiT[v^] / {Skg) and r = 300 
K the input equilibrium temperature. Results correspond to 
three-dimensional simulations of water at p = 0.6 g/mol/A^ 
with spatial resolution 5 = 20 A in each direction. 



The accuracy of the stochastic time integrator for 
the Langevin equation affects the value of the numer- 
ical temperature. Figure [5] shows the dependence of 
the relative error in the mean temperature (defined as 
Et = \Tnurn - Te\/Te) with the time step At. Good 
agreement is found and the relative error remains smaller 
than 10% for At < 100 fs. 

Figure [6] shows the standard deviation of the cell mass 
density Std[pc] against the mean density p — M/V, 
where M and V are the total mass and volume of the 
system. The grand canonical prediction for the equi- 
librium state is Std[p] = [pksT / {c^^Vc)]'^/'^ and IS com- 
pared with the numerical simulations. In the ideal gas 
limit the sound velocity is just ^JksT /m (with m the 
molecular mass) and the density fluctuations increase as 
Var[p^*'''^°')] = {m/V)p. As the fluid becomes denser, it 
becomes less compressible (the isothermal sound speed 
ct increases) and as a consequence in the liquid phase 
the mass fluctuations decrease substantially. Therefore, 
the largest mass fluctuations are observed at moderate 
densities (e.g. around p ~ 0.3 g/m.o\/K? for argon 

see Figl6]) . Almost perfect agreement is found between 
theoretical and numerical results for both argon and wa- 
ter. 



B. Correlations at equilibrium 

We now show that fluctuations are transported 
through the system and dissipated in a hydrodynamically 
correct way. This can be shown via the time correlation 
of the fluctuating quantities. As shown in Ref. the 
time correlation of the Fourier components of mass and 
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FIG. 6: Standard deviation of the mass density S[p\ in a fluid 
cell with volume 37.5 nm'^ of argon and water at T = 30QK. 
The continuous line corresponds to the grand canonical results 
obtained using the equation of state for argon, the dotted lines 
for water (TIP3P model) [pkeT / {c%Vc)Y^'^ , and the dashed 
lines show the ideal gas limits [(m/I/)p]^/^ with niAr = 39.498 
g/mol (upper curve) and niH^o = 18.015 g/mol (lower curve). 
Circles are results from the fluctuating hydrodynamics solver 
(using a water simulation at p = 0.632 g/mol/A"^). 



momentum satisfy 

(p(fc,t)p(fc,0)) 
Var[p(A:,0)] 
(g||(fc,t)ff||(fc,0)) 

Var[5||(fc,0)] 
(gx(fc,t)ffj.(fc,0)) 

Var[gi(fc,0)] 
(p(fc,^)^g||(fc,0)) 
(p(fc,0)zg||(fc,0)) 



exp{— Frfc^t} cos{cTkt), 
cxp{— Frfc^i} cos(cTfct), 
exp{— i/fc^t}, 

exp{-FTfc^t}sin(cTfci). (14) 



In Eqs. p4)) indicates the longitudinal momentum, par- 
allel to the wave vector fc, while g± indicates the transver- 
sal components. Sound modes couple longitudinal mo- 
mentum and density, as shown in Eas. fTi]) . According 
to the Landau description of fluctuating hydrodynam- 
ics PI, at a fixed time the stress fluctuations occurring 
within different fluid volumes are spatially uncorrelated. 
This means that the variance of the Fourier modes of 
the hydrodynamic variables (Var[p(/c, 0)] = (p(/e,0)^)) is 
independent of the wave-number fc. Moreover, fluctua- 
tions from the equilibrium state are assumed to be small 
so that a linear analysis can be applied. This means 
that perturbations of different wavelengths evolve inde- 
pendently as correlations between fluctuations with dif- 
ferent wave- vectors are negligible. 

We have evaluated the time correlations of the Fourier 
components of the hydrodynamic variables in a periodic 
domain for the set of wave- numbers fc„ = 2i:n/L. These 
correlations were then fitted to the corresponding expo- 
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FIG. 7: Water at p = 0.632 g/mol/A^ and T = 300 K (am- 
bient conditions) within a periodic box of size 50 x 50 x 2250 
at equilibrium; the FH mesh is comprised of of 1 x 1 x 150 
cells, (a) The density field in the real space, p[x,to). (b) 
The time-dependence of the density at one cell p{xo,t). (c) 
The Fourier mode p{ki,t) associated with wavevector fei = 
{0,0,2n/Lz). (d) The (normalized) time correlation function 
(p(fei, i)p(fci, 0)). Note the different time scales associated 
with variations of quantities in (b),(c) and (d). The entire 
run is of 40 ns duration. 



nentially decaying functions of Eqs. p^ to obtain the ef- 
fective sound frequency crk and the effective decay rates 
for each wave-number (i.e. the inverse of Trk^ and vk^ 
for sound and shear, respectively). Calculations were 
done in a periodic system of size 50 x 50 x 2250 A'^, 
with a mesh of 1 x 1 x 150 cells, and we considered per- 
turbations with wave-vectors fc„ = (0,0,27™/^^). In or- 
der to illustrate the mathematical transformations, Fig[7] 
shows the density field in the real space p(x,to) (a), the 
time-dependence of the density at one cell p{xo,t) (b), 
the Fourier mode p{ki,t) associated with fci = 27r/i^ 
(c) and, finally in (d), the time correlation function 
{p{ki,t)p{ki,0)) together with the best fit obtained to 
the theoretical exponential decay of the sound mode. The 
best fits to the effective decay rates and sound periods 
for varying wavelength A„ = 27r/fc„ are compared with 
the theoretical relations in FiglH Note that the theoreti- 
cal trend agrees quite well with the simulation results for 
A > 100 A. Considering that the size of one cell in these 
simulations is 15 A, this means that one needs more that 
about 7 cells to properly resolve one wave. In fact, for 
wavelengths A < 100 A the viscosity is underestimated 
due to the reduction in spatial resolution. The same rea- 
soning applies to the sound time. 



V. NON-EQUILIBRIUM STATES 

In this section we present standard non equilibrium 
flow tests performed with the mesoscopic Eulerian solver 
described in previous sections. In order to do so, firstly 
we require the description of the explicit boundary con- 
ditions for our fluid system enclosed between walls. 



A. Boundary conditions 

The imposition of boundary conditions on the velocity 
is illustrated in FigEl In the figure the boundary "li;" 
is placed at the interface between cells and 1. The 
fluid region corresponds to cell and continues to the 
right, while cells 1 and 2 (in grey) are ghost cells which 
are used to impose the desired mechanical behaviour at 
the boundary "w" . In order to close the system, we 
need to evaluate the momentum flux across the inter- 
face w. As at any other cell interface, we approximate 
n.u, • = [(IIo 4- ni)/2] • Cu,, where e^, is the surface 
unit vector (in this case eoi). Hence we require knowl- 
edge of III, the stress tensor in the first ghost cell as 
can also be inferred from Eq. ([6]) when evaluated at the 
boundary cell "0" . According to the constitutive relation 
Hi oc Vvi. Hence, a secondary ghost cell (#2 in FigE]) 
is required to evaluate the velocity gradient via the cen- 
tral difference scheme: Vui • e^, = (2S)~^{vq — V2) ■ e^, 
(where S is the spatial resolution). This closes the set of 
equations for the velocity at the ghost wall cells. This pro- 
cedure enables certain flexibility: one can either impose 
the value of the momentum flux 11^, • (von Neumann 
boundary condition) , a generalized relation involving the 
fluid velocity at the wall 11^ • e^, oc v^, (Maxwell relation 
for fluid slip [37*1) or the more standard no-slip condi- 
tion, v-uj = Uwaii (Dirichlct boundary condition). In the 
present work we assume no-slip at the wall and set the 
velocity of the ghost cells accordingly to a linear interpo- 
lation of the velocity; with ~ U^aii we find 



^^1 

V2 



- 2 Uwall — ■UQ, 



(15) 



More general slip boundary conditions can be obtained 
by choosing a different value for As is customary, 
the density at the wall is uniquely controlled by the fluid, 
meaning that M„ = Mq = Mi = M2 [2^ . 



B. Couette, Poiseuille and cavity flows 

Stochastic and deterministic simulations for three dif- 
ferent flow situations (Couette, Poiseuille and cavity 
flow) have also been performed, displaying good com- 
parisons with theoretical predictions. Figures [TOl and [TT] 
show stationary fluctuating (and deterministic) flows for 
argon at ambient temperature and mass density po = 
0.599 g/mol/A'^. The simulations are performed using 
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FIG. 8: The best fits to the decay times obtained from the time correlations of the Fourier modes of wavelength A — 2n/k. 
Symbols correspond to numerical results and lines to theoretical predictions (see Eqs. (|14p '). (a) Shear decay times: the lines 
are the theoretical results {k^v)~^, (b) sound attenuation time (/c^Ft)"^ and (c) the sound period \/ct- Circles are obtained 
from the imaginary part of {h{k,T)h{k,0)) and squares from the real part, where h — v± (transversal velocity) in (a) while 
ft = p in (b) and (c) (similar results, not shown, were obtained for the longitudinal velocity iiy). Results for water (ct = 14.75 
A/ps, u — 84.98 A'^/ps and C/p ~ 201.02 A^/ps) correspond to the simulations illustrated in Fig[7l while argon {ct ~ 5.61 
A/ps) corresponds to p = 0.6 g/mol/A'^ and T = 300 K in the same box used in Fig[7] 
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FIG. 9: Boundary conditions imposed at the interface "w" . 
Cell represents a fluid cell, and cells 1 and 2 are ghost cells. 
Arrows indicate the velocity vectors in a case with zero ve- 
locity imposed at "w" (no-slip boundary condition and rigid 
wall at rest). 



10x10x10 cells representing a periodic box of size 
200 X 200 X 200 A^. The fluid is confined in a chan- 
nel defined by two infinite parallel planes orthogonal to 
the z axis. In this particular geometry the first two and 
last two layers of particles in the z direction are ghost 
particles belonging to the walls; the no-slip condition is 
satisfied at = 30 A and Z2 — 150 A, so the fluid is 
confined in a region of width 120 A. 

A Couette flow is shown in Fig[Tni We plot the x 
component of the stationary velocity field in the z direc- 
tion. The wall amplitude velocity has been set at 2.04 
A/ps, while the amplitude of fluid velocity fluctuations 
is about 0.5 A/ps, for the temperature and cell volume 
considered. The inset picture corresponds to an equiva- 
lent simulation but with the thermal fluctuations within 
the pressure tensor switched off. In this limit, we recover 
standard Navier-Stokes behaviour. 

Figure [TT] shows the x component of the stationary 
velocity field in the z direction for a Poiseuille flow. The 
applied gravity force in this case is 0.0174 A/ps^. The 
same simulation performed without thermal fluctuations 
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z(A) 

FIG. 10: Stationary Couette profile according to fluctuating 
hydrodynamics: The diamonds are simulation results and the 
continuous line is the theoretical stationary linear profile. The 
wall velocity has been set 2.04 A/ps. The inset figure shows a 
deterministic Navier-Stokes simulation with the same param- 
eters. Vertical dashed lines represent the boundary walls. 



is displayed in the inset picture and both are compared 
with the theoretical solution. 

We have also carried out simulations of cavity flow for 
(TIP3P) water depicted in FiglT^ within a domain of di- 
mensions 1500 X 1500 X 50 A^ and mesh 30 x 30 x 1. The 
wall is moving at a constant speed of 1 A/ps in the y di- 
rection. The average fiow corresponding to averaging the 
fluctuating hydrodynamics result (Fig. 112b ) corresponds 
to the fluctuation-less fiow of Fig. [T2b . Although this is 
quite a large fluid domain, thermodynamic fluctuations 
arc still very visible and have a major effect on the flow. 
For smaller domains the stationary circulatory fiow can 
be completely nullified by the fluctuations. This kind 
of cavities can be used as mixers in microfluidics appli- 
cations, but the extent of the mixing is affected by the 
magnitude of the ffuctuations. 
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tuations which are important at the mesoscale, typically 
in the nanometer range, as well as the hydrodynam- 
ics. We have tested the model in equilibrium and non- 
equilibrium situations. Simple no-slip boundary condi- 
tions have been used and tested in three flow situations 
Couette, Poiseuille and cavity flow. From a set of molec- 
ular dynamics simulations, we have also derived a simple 
approximation to the equations of state for the TIP3P 
water model and for argon which permit us to simulate 
these compressible fluids at ambient temperatures ( ca 300 
K) and pressures around 1 atm. 



FIG. 11: Stationary Poiseuille profile according to fluctuating 
hydrodynamics: Diamonds correspond to simulation results 
and the continuous line is the theoretical profile associated 
with Navier-Stokes fiow. The inset figure shows the same fiow 
case for a purely deterministic simulation. In both figures, the 
vertical dashed lines represent the walls. 
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There are some aspects of the model which could be 
improved and we reserve for future work. The discretiza- 
tion of the second order derivatives is based on a central 
difference method. This approach is known to produce 
instabilities which lead to a characteristic and undesir- 
able alternating pattern in the velocity and density fields 
[38| . However, we have found that such problems only 
arise with very strong perturbations well beyond the typ- 
ical flows at these scales. A further important extension 
is to incorporate the energy equation so as to be able 
to model thermal phenomena. Finally, we have used a 
very simple boundary condition appropriate for the cases 
studied here: more sophisticated boundary conditions, 
which may reproduce diverse molecular boundaries (e.g. 
taking into account the hydrophobic hydrophilic nature 
of the walls) while also retaining the fluctuations, can be 
devised in the framework of the hybrid models [j] . 

Our mesoscopic fluctuating model serves to support a 
wide range of applications. We are currently using it 
in hybrid molecular-continuum simulations [4], with im- 
plicit solvent models Q, and we plan to use it for the 
study of microfluidic flows. The use of a regular lattice 
and the Eulerian description greatly simplify the imple- 
mentation of this model in a serial (as here) or parallel 
computing environment. As a result, we believe that it 
furnishes a unique tool to explore hydrodynamics at the 
nanoscale including the effects of fluctuations, in stand- 
alone more or coupled with molecular dynamics. 



FIG. 12: Streamline plots of a cavity flow for TIP3P water 
model in the stationary regime. The wall velocity along y is 
1 A/ps, and the temperature is 300 K. The dimensions of the 
cavity are 1500 x 1500 x 50 and the mesh is 30 x 30 x 1. 
In (a) thermal fluctuations of the pressure tensors are not 
present, and the Navier-Stokes solution recovered. In (b) we 
show the fluctuating hydrodynamics solution. 



VI. SUMMARY 

We have derived a finite volume discretization of the 
equations of fiuctuating hydrodynamics. The model pro- 
vides a good representation of the thermodynamic fiuc- 
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APPENDIX A: HYDRODYNAMIC MODES 

In the isothermal situations the equations that describe 
the transport of mass and momentum density fields are 

dtp = -V • g, 

^tg = -V • (gv) - Vp + + (I + c) V (V • v) . 

(Al) 

The equilibrium state is characterized by a constant den- 
sity field Pe and a zero momentum density field gg = 
because the fluid is at rest. For fluctuations of small 
enough amplitude, the relaxation towards equilibrium is 
governed by the linearized version of the mass and mo- 
mentum Eas. (|Aip . By decomposing the hydrodynamic 
fields as p(r, i) = + 5p{v, t) and g(r, t) = 5g(r, i), the 
linearized version of Egs. ljAip for the perturbations re- 
sults 

dtdp = -V • (5g, 

dt5^ = -Vp + vVHg + UBV {V -Sg), (A2) 

with the usual definitions for the kinematic shear vis- 
cosity v = rj/pe and the effective bulk viscosity vb = 

(77/3-f C)/Pe. 

Under isothermal conditions the thermodynamic rela- 
tion of pressure perturbation Sp with density and tem- 
perature becomes simply [39| . 

6p{p,T) = cUp)5p, (A3) 

where = {dp/dp)T is the squared isothermal sound 
velocity. 

Let us consider a general solution as a series of normal 
modes 

a(r,t) = a(fc,^)e*'=■^ 

where we have gathered the hydrodynamic variables in 
the array a — {Sp, Sg}. Taking the Fourier transform in 
Eqs. ljA2[) one gets the equations a{k,t). In the linear 
regime there is no coupling between modes and without 
loss of generality one can work in the reference system 
for which the wave vector is fe = (fc, 0, 0), 



da{k, t) 
Jt 



= Ha{k, t) 



where a{k,t) = {p{k,t), gx{k,t), gy{k,t), g^{k,t)) and 
the hydrodynamic matrix is 



H 



ik 

ic^k i^Lk^ 



lyP 



(A5) 



with the kinematic longitudinal viscosity defined as = 
p + pb- The eigenvalues of the hydrodynamic matrix 



H provide the growth rates of the normal modes of the 
system given by Eq. (|A4p . The eigenvalues are obtained 
from the roots of the characteristic equation det[H — 
Lul] ~ 0, which results in 

(w + z^fc2)^(w2-f^i/Lfc^w + 4fc2) =0. (A6) 
The solutions are 



wi,2 = -vk , 

= —Yxk^ ± isTk, 



(A7) 



where we have defined as the isothermal sound ab- 
sorption coefficient and st as the sound speed depending 
on the wave vector given by 



Ft = 



PL_ 

2 



St = 



(A8) 



The first two eigenvalues (^1,2) correspond to the two 
shear modes associated with the exponential decay of the 
transversal momentum gy and gz- Sound modes corre- 
spond to a;3^4. Indeed, as can be seen from Eqs. (jATj 
and Eqs. (|A8|) . sound is underdamped if st is a real 
number. However, according to Eq. (|A8p ii k > 2ct/vli 
sound becomes overdamped. Nevertheless for most liq- 
uids ct/pl ^ 0{1) so this anomalous solution occurs at 
quite small wavelengths for which the present mesoscopic 
description does not apply (at molecular lengthscales one 
should consider the dependence of the transport coeffi- 
cients on k within the generalized hydrodynamic formal- 
ism [3^). As a matter of fact, the difference between st 
and ct is negligible for any mesoscopic wavelength, so 
throughout the present paper we assume that st = ct- 
With this last approximation, the solution is given by 

p{k,t) — pe + p{k,0) expl—Trk^t} cos{cTkt) 

- — exp{-FTfc^t}sin(cTfci)A;-g(k,0), 

CT 

g{k,t) = exp{-FTfc^O[cos(cTfci)fc • g(A;,0) 

— ism{cTkt) CTp{k,0)]k 

+ exp{-iyk'^t} (1 - kk) ■ g(fc, 0) (A9) 



(A4) where k = k/\k\ is the unit wave vector. 



APPENDIX B: EQUATIONS OF STATE FOR 
ARGON AND WATER VIA MOLECULAR 
DYNAMICS SIMULATION 

In this section we study the equations of state for water 
and also argon through molecular dynamic simulations 
using the NAMD molecular dynamics code [1^. In par- 
ticular, the theoretical Lennard- Jones equation of state 
for argon given in Refs.flO, [Ml, [s^l is not necessarily ex- 
act because our MD simulations using the CHARMM 
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FIG. 13; Argon equation of state: Pressure (in bars) versus 
mass density at temperature 300K. The simulation results 
appear with error bars and the dashed line is the equation 
of state for the theoretical Lennard- Jones fluid in Ref . [30l. [3ll. 
[3^ . The continuous line is the fit of the numerical data to 
the second order polynomial 3088.21 - 12065. 2p + 14765. 8p^ 
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FIG. 14: Water pressure (in bars) versus temperature ob- 
tained from MD simulations of the TIP3P water model (using 
NAMD) at a fixed density p = 0.55066 g/(mol A^). The * 
symbols correspond to simulations with all bonds rigid while 
the □ symbols are for non rigid bonds. 



force field perturb the Lennard- Jones potential close to 
the cutoff radius, smoothing it to zero |40|]. Thus we de- 
cided to obtain an accurate approximation of the equa- 
tion of state by directly fitting the data obtained from 
molecular dynamics simulations to a second order poly- 
nomial. These considerations also apply for the TIP3P 
water model. 

We have computed the equation of state for argon in 
our molecular model in the NpT ensemble at a fixed tem- 
perature 300 ± 4 K with 5000 argon atoms. The simula- 
tions are performed with a time step At — 1 is. We have 
made use of a switched Lennard-Jones potential [40] be- 
tween 10 and 12 A. The results are displayed in Fig [13] 
for the mass density. Note that the deviations between 
the theoretical model and the simulation results for the 
pressure can be as large as 15% in the units presented in 
the graph, the simulated argon pressure always being less 
than that given by the theoretical Lennard-Jones equa- 
tion of state in Ref. [3^, i31„ ^] . 



Concerning the water equation of state, we have per- 
formed simulations of water molecules using the TIP3P 
water model |2|]. In Fig[T3]we present the simulation re- 
sults for a range of temperatures around 300K with a 
time step At = 1 fs. We have chosen a fixed density 
p — 0.55066 g/(mol A^) what gives (in a NVT ensemble 
for a cubic periodic box of size 30 A) a total of 826 water 
molecules. We have also tested the effect of using rigid 
bonds. We see that for harmonic bonds the fluctuations 
in pressure are bigger compared to those of the rigid sim- 
ulation providing also smaller mean average pressures in 
general terms. We have also compared our results with 
three analytical models for the water equation of state 
presented in Ref. [s^. We observe clear deviations from 
the models. Basically, all three theoretical models largely 
overestimate the isothermal sound velocity at any value 
of the density considered. This reinforces the necessity of 
pre-calibrating the equation of state for each particular 
fluid considered; for instance via MD simulations as done 
here. 



In Fig[T5]we plot the values of mass densities obtained 
from water MD simulations in an NpT ensemble against 
pressure (in bars) at a fixed temperature 300 ± 5 K with 
non-rigid water molecules and a time step At = 1 fs. The 
second order polynomial fit used in the FH equations is 
also shown in the continuous line. 
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